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Over the last century, a large number of physical and mathematical developments paired 
with rapidly advancing technology have allowed the field of quantum chemistry to ad- 
vance dramatically. However, the lack of computationally efficient methods for the exact 
simulation of quantum systems on classical computers presents a limitation of current 
computational approaches. We report, in detail, how a set of pre-computed molecular 
integrals can be used to explicitly create a quantum circuit, i.e. a sequence of elementary 
quantum operations, that, when run on a quantum computer, to obtain the energy of 
a molecular system with fixed nuclear geometry using the quantum phase estimation 
algorithm. We extend several known results related to this idea and discuss the adia- 
batic state preparation procedure for preparing the input states used in the algorithm. 
With current and near future quantum devices in mind, we provide a complete example 
using the hydrogen molecule, of how a chemical Hamiltonian can be simulated using a 
quantum computer. 

Keywords: electronic structure, quantum computing 



1. Introduction 

Theoretical and computational chemistry involves solving the equations of 
motion that govern quantum systems by analytical and numerical methods [H 
[2]. Except in standard cases such as, the harmonic oscillator or the hydrogen 
atom, analytic solutions are not known and computational methods have been 
developed. 

Although classical computers have tremendously aided our understanding of 
chemical systems and their processes, the computational cost of the numerical 
methods for solving Schrodinger's equation grows rapidly with increases in 
the quality of the description. Research is ongoing to improve computational 
methods, but large molecules and large basis sets have remained a consistent 
problem despite the exponential growth of computational power of classical 
computers [3j. 

Theoretical computer science suggests that these limitations are not mere 
shortcomings of the programmers but could stem from the inherent difficultly 
of simulating quantum systems. Extensions of computer science using quantum 
mechanics led to the exploitation of the novel effects of quantum mechanics 
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for computational purposes resulting in several proposals for quantum com- 
puters Quantum simulation is the idea of using quantum computational 
devices for more efficient simulation (SJ IH] • Since the dynamics are simulated 
by a quantum system rather than calculated by a classical system, quantum 
simulation often offers exponential advantage over classical simulation for cal- 
culation of electronic energies [7] , reaction rates [HI E] , correlation functions [TU] 
and molecular properties . Recently, a review of these techniques and other 
applications of quantum computing to chemistry has appeared [I2] . 




This paper does not consider the effect of errors, however it is an important 
consideration that needs to be taken into account. Quantum error correction 
methods have been developed to counteract the unwanted effect of quantum 
noise, however fault tolerant constructions require redundant qubits and only 
allow a discrete set of gates to be used [131 E] • This is not a serious cause for 
concern as the conversion from a continuous set of gates to a discrete set of 
gates only requires a poly-logarithmic overhead [15] . Clark et al. [T6] estimated 
the resources required to compute the ground state of a one dimensional trans- 
verse Ising model and found, using experimental parameters from a proposed 
ion trap quantum computing implementation, that the fault tolerant construc- 
tions would be too costly for straight-forward applications to simulation. This 
suggests that quantum simulation without quantum error correction is more 
feasible for the near future. 

The state of the art in experimental realizations of quantum simulation for 
chemistry is represented by calculations of the energy spectrum of molecular 
hydrogen first by Lanyon et al. [17] using an optical quantum computer. Very 
soon after, Du et al. [IB] used NMR technology to demonstrate the adiabatic 
state preparation procedure suggested by Ref. [^ as well as reproduce the 
energy to higher accuracy. 

The key limitation of both experimental algorithms was the representation 
of the simulated system's propagator. Both experiments, relied on the low 
dimensionality of the propagator for the minimal basis H2 model considered. 
The unitary propagator for a two-level system can be decomposed using the 
real angles a, /3, 7: 



Due to this decomposition, longer propagation times corresponding to higher 
powers of unitary evolution operator, can be achieved by changing a to 
ja and 7 to 77, thereby avoiding the need of further decomposition of the 
unitary operator. Beyond the two dimensional case, this decomposition is not 
available. 

The objective of this paper is to provide a general decomposition for elec- 
tronic Hamiltonians and demonstrate this method with an explicit quantum 
circuit for a single Trotter time step of the minimal basis hydrogen molecule. 
This is an extension of the supplementary material from Lanyon et al. |17| . 

2. Overview of the quantum algorithm 

The use of the Fourier transform of correlation functions in computational 
chemistry allows for the extraction of information about many properties such 
as transport coefficients |19j and molecular spectra [20] . For the application of 
quantum computation to molecular electronic energy calculations, the spec- 
tral solution of the Schrodinger equation [21] is pursued. The key idea of the 
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Figure 1. An algorithmic overview of the steps taken to simulate a chemical Hamiltonian on a 
quantum computer. The time independent Hamiltonian of a molecular system (to the left of the box) 
will be decomposed into a sum of Hermitian matrices, (hi) and by means of a Trotter decomposition 
the unitary propagator U can be constructed. The Jordan-Wigner mapping is used to convert the 
propagator into a sequence of quantum gates. Phase estimation algorithms are used to evaluate the 
eigenvalue of a correctly prepared stationary state \ipo). 



approach is that the Fourier transform, J^[f{t)] = J exp{iujt)f(t)dt of the 
autocorrelation of the time evolving state, P{t) = (^(0)|^(t)), has peaks at 
the eigenenergies. Typically, the semiclassical propagator is used for examin- 
ing vibrational structure e.g. |22j . Since quantum computers can simulate the 
quantum evolution of electrons efficiently, the Fourier transform of the elec- 
tronic autocorrelation function can be calculated directly. The measurement 
technique presented here collapses an input wave function into an eigenstate 
and returns the frequency of the autocorrelation of the time-evolving eigen- 
state by employing the spectral method mentioned before to the extraction 
of molecular electronic energies. In this paper, we focus on the extraction of 
eigenvalues given input states using the time evolution of the state. 

The construction of the general quantum circuit to simulate the evolution 
of the molecular system is performed in three steps: 

(1) Write the Hamiltonian as a sum over products of Pauli spin operators 
acting on different qubits. This is described in Section [3] and made possible 
by the Jordan-Wigner transformation. 

(2) Convert each of the operators defined in step (1) into unitary gates 
such that their sequential execution on a quantum computer can be made 
to recover an approximation to the unitary propagator, exp{—iHt). This 
is detailed in Section HI 

(3) The phase estimation algorithm, as described in Section [5| is then used 
to approximate the eigenvalue of an input eigenstate using the quantum 
Fourier transform of the time domain propagation. Section [6] discusses 
eigenstate preparation. 

To demonstrate these steps, the construction is applied to the example of 
the hydrogen molecule in Section [7j The key components of the simulation 
procedure are depicted in Fig. 1. The next section provides a basic review 
of the fundamental concepts and notations of molecular quantum chemistry 
for the benefit of quantum information scientists and to establish the nota- 
tion. A detailed account of electronic structure methods can be found in the 
monographs [TJ |2] . 



3. The Electronic Hamiltonian 

The Born-Oppenheimer approximation assumes that the wave function of a 
molecule can be expressed as a product of the nuclear wave function and 
electronic wave function (parameterized by the nuclear coordinates) due to 
the difference in mass between electrons and nuclei. This approximation allows 
for the solution of the time independent Schrodinger equation of the electronic 
wave function for a given nuclear geometry. 
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Table 1. An overview of second quantization for Fermionie particles. For quantum chemistry, the annihilation and creation oper- 
ators correspond to removing (adding) an electron into a particular orbital, {xfc}- The anti-symmetry is enforced by the canonical 
commutation relations, and the A^-clectron wave function is expanded over the configuration state functions of the Fock space. 



Second quantization 

Creation operator a||ji, . . . , 0,;, . . . , j„) = r^lji, 1^, j„) with = n^=\(-l)^ 

a\\ji,. . . , li, . . . ,jn) = 
Annihilation operator ai\ji, . . . . . . , j„) = P||ji, 0^, Jn> 

ai\ji,...,Oi,...,jn) = 
Canonical commutation relations {aj,aj} = Sijl 

{ai.aj} = 



Fock space 
Basis vectors, 

configuration state functions 

Inner product 
Vacuum state 

Single electron operator, A{x) 



■■■Jn) 

n 



p=i 

(j|k> = nf=i5ip,fcp 

{vac\vac) = 1 
a; \vac) = 



\vac) 



where ji = 0, 1 



where Aij = f Xi(xi)M^l)Xj(^l)d^l 
is a single electron operator 
and {xk} corresponds to {a^} 



The molecular electronic Hamiltoniarj^ in second quantized form is given by 

mm- 

H = ^ ^ hpqCLpCLq + — ^ ^ hpqrsOipCL^gClrO'Sj (1) 
p,q p,q,r,s 

where the sum is over the single particle basis set described below. The annihi- 
lation {ttj} and creation operators {a]} obey the Fermionic anti-commutation 
relations (see Table [l]): 

[ap, ag]+ = 0, [ap, aj]+ = 6pgl, (2) 

where the anti-commutator of operators A and B is defined as [A, = 
AB + BA and 1 is the identity operator. The annihilation (creation) operators 
correspond to a set of orbitals, {xi}-, where each orbital is a single-particle 
wave function composed of a spin and a spatial function, denoted cTj and 
respectively. As the Hamiltonian considered here commutes with the electron 
spin operators, (Ti is restricted to be one of two orthogonal functions of a spin 
variable uj that we denote a{uj) and /3(u;). Similar Hamiltonians can be found 
in many physics problems involving Fermionic particles. 

Although any basis can be used, the molecular orbitals are particularly con- 
venient for state preparation reasons discussed below. The molecular orbitals, 
in turn, are formed as a linear combinations of atomic basis functions |23|l24j. 
The coefficients of this expansion are obtained by solving the set of Hartree- 
Fock equations which arise from the variational minimization of the energy 
using a single determinant wave function. Due to its restriction to a single 
determinant, the Hartree-Fock solution is a mean field solution and the dif- 
ference between the Hartree-Fock solution using an infinite basis of atomic 



^Throughout this paper, atomic units are used: h (1.054 -10 '^^ J s), the mass of the electron 
(9.109 ■ 10"'^^ kg), the elementary charge (1.602 ■ 10~^^ C), and the electrostatic force constant 
(l/47reo = 8.988 ■ 10^ N m^ C'^) are set to unity. 
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orbitals and the exact (correlated) solution defines the electron correlation 
energy. 

The matrix elements {hpq\ and {/ipgrs} hi Eq. ([T]) denote the set of one- and 
two-electron integrals that must be evaluated using a known set of basis func- 
tions (the basis set) during the Hartree-Fock procedure. Ideally, the number of 
basis functions used would be infinite, but in practice only a finite number of 
basis functions are used. By selecting Gaussian single-particle basis functions, 
these integrals are efficiently computable. Next, to further establish notation, 
we develop the explicit form of the integrals hpq and hpqrs- 

We denote the set of single-particle spatial functions which constitute the 
molecular orbitals {'^k{x)}^^^. Finally, define the set of spin orbitals as 
{Xp{'^)}lp=i with Xp = Vi'^i aiid X = (x,w) where cjj is a spin function. The 
one-electron integrals involving the electronic kinetic energy and the electron- 
nuclear attraction terms are: 

V ^ / dxx;(x) ^-iv2 - ^ ) x,(x) (3) 

and the two-electron integrals involving the electron-electron interaction, 1 /ri2 
are: 

U [a A X;(xi)Xg(x2)Xr(x2)Xs(xi) 

hpqrs = / dxidX2^^^ (4) 

J n,2 

In Eq. ([3]), is the Laplacian with respect to the electronic spatial coor- 
dinates. The positive valued scalars r^^x and ri 2 are the Euclidean distance 
between the a**^ nucleus and the electron and the Euclidean distance between 
the two electrons. In both Eqs. ^ and Q, the spin of the electron can be in- 
tegrated out resulting in integrals over the spatial components however since 
the focus is on spin orbitals, the electron spin will be included in the definition 
of the integrals. 




3.1. Representing the molecular Hamiltonian in terms of quantum bits 



Just as classical computation is based on the notion of a bit, the basic unit 
of quantum information is the quantum bit (qubit). Qubits, being quantum, 
are described by a wave function instead of a probability distribution as in the 
case of classical bits. The use of the wave function description allows for the 
superposition of states and for entanglement between different qubits. More- 
over, since qubits themselves are described using wave functions, quantum 
states can be efficiently stored. 

In principle, any two-level quantum mechanical system can be considered 
a qubit. Examples include photons, ions, and super conducting loops. Practi- 
cal requirements for qubits and their manipulation was originally outlined by 
DiVincenzo [25] and experimental progress towards satisfying the DiVincenzo 
criteria were recently reviewed [3]. Since two- level systems can be described 
as spin-half particles, the relevant (Pauli) matrices are: 



a 



1 

1 



-i 

1 



a 



1 
-1 
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Notice that we have defined the Pauh matrices with eigenvalues ±1 instead 
of Together with the identify matrix, 1, the Pauh matrices form an 

operator basis for two-level systems. The eigenvectors of cj^ are labeled as 
|0) and |1) corresponding to the eigenvalues +1 and -1 respectively. There are 
several computationally equivalent models of describing quantum computation 
but here we focus on the circuit model of quantum computation. A more 
comprehensive introduction to quantum computation can be found in Ref. [13] . 

The quantum circuit model uses a set of elementary gates to reproduce the 
action of arbitrary unitary transforms. If the elementary gates are capable of 
reproducing any desired unitary transform to arbitrary precision, the set of 
gates is called universal. A universal gate set requires single qubit gates and 
any two-qubit entangling gate. The two-qubit CNOT gate leaves one qubit 
space unchanged and acts with on the second qubit when the first qubit is 
in the state The gate is given as CNOT= |1)(1| ® cr^ + |0)(0| (g)l. The set, 
Rx, Ry, and Rz gates can generate any single qubit gate where Rn is defined 
as exp[— i(j"^/2] for real 9. 

For experimental addressability, the qubits must, in general, be distinguish- 
able. However, the electrons of the molecular system are indistinguishable. 
The Jordan- Wigner transform is used to circumvent this issue by expressing 
Fermionic operators in terms of the Pauli spin operators {a^, o"^, o"^, 1} that 
correspond to the algebra of distinguishable spin 1/2 particles [101126]. The 
Jordan- Wigner transform is given by: 





"l o" 
1 


®i-i 


"O l" 






1 




o" 
-1 


^N-j- 




"l o" 
1 


Cg)j-1 


"o o" 

1 




1 




o" 

-1 





1 

(5a) 



where a+ = ^^^±^ = |0)(1| and a' = ^^^^ = |1)(0|. The qubit state |0...0) 
corresponds to the vacuum state and the string of operators, preserve the 
commutation relations in Eq. ^ since and anti-commute. The spin 
variable representation of relevant operators after the Jordan- Wigner trans- 
formation is given in Table |X2j 

4. Efficient approximations of the unitary propagator by a Trotter 
decomposition 

As mentioned in the introduction, the idea of using a quantum computer 
for obtaining molecular energies requires efficiently approximating the uni- 
tary propagator, exp{—iHt), for a time sufficiently long to resolve the Fourier 
frequency to a desired precision. The present section continues by describing 
the Trotter decomposition for non-commuting Hamiltonian terms and then 
presents the quantum circuits for the relevant exponentials. 

4.1. Trotter decomposition 

Using the second-quantized representation allows for a straightforward decom- 
position of the exponential map of each term of the Hamiltonian. However, the 
terms of this decomposition do not always commute. The goal of the Trotter 
decomposition is to approximate the time evolution operator of a set of non- 
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commuting operators. The operators are exponentiated individually for small 
time steps and the procedure is repeated such that their product provides a 
reasonable approximation to the exponentiation of the sum. Using this approx- 
imation, the construction of the time propagator can be efficiently carried out 
on a quantum computer provided that the Hamiltonian can be decomposed 
into a sum of local Hamiltonians [B] . The first-order Trotter decomposition is 
given by: 



-iHt 



t 



where t/At is the Trotter number [27]. As the Trotter number tends to infinity, 
or equivalently At — )• 0, the error in the approximation vanishes. If t/At is not 
an integer, the remainder is simulated as another Trotter time slice. There exist 
higher order approximates (Suzuki- Trotter formulas) which reduce the error 
of approximation even further. For instance, the second order approximations 
is given by: 



^) e-/.«At . . . e-^'^^^)) ^'+0{t{Atf 

(7) 

Higher order approximations take increasingly more complicated forms |27] 
and were first studied in the context of black box quantum simulation of sparse 
Hamiltonians by Berry et al. [S^. They considered Hamiltonians composed of 
m efficiently implementable terms and showed that the number of exponen- 
tials cannot scale better than linear in the time desired and the maximum 
frequency of the full Hamiltonian. The proof shows that if sublinear simula- 
tion of arbitrary Hamiltonians were possible, bounds for the power of quantum 
computation proven in Ref. [29] could be violated leading to a contradiction. 



4.2. Quantum circuit primitives 

Each term of Eq. ([T]) can be exponentiated using the universal gate set de- 
scribed in Subsection |3.1| after performing the Jordan- Wigner transformation 
to Pauli spin matrices. We will outline the procedure for generating quantum 



circuits for chemical systems and summarise the results in Table Al The con- 
struction of quantum circuits for general Fermionic Hamiltonians is further 
discussed in Ref. [I0l[30l[3l]. 

To understand the exponential map of the product of Pauli spin matrices, 
first consider the exponential map of two operators. To create the unitary 
gate exp[— i|((T^(8)(T^)], the CNOT gate can be used to first entangle two qubits, 
then the Rz gate is applied, and followed by a second CNOT gate fT3]. 



-< 




Rz 


-< 













This construction can be generalized to more qubits by using additional CNOT 
gates. For example, the circuit for the three-body operator involving three 
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qubits, exp[— i|(o-^ ® o-^)], is simulated by the following quantum circuit: 



1 


1 


< 

k A r 


>— 

y- 




r 




Rz 

















As seen from the three-qubit example above, this construction can be readily 
extended for n-fold products of operators. 

4-2.1. Construction of different Pauli matrix tensor products 
If one requires a different tensor product of Pauli matrices besides the prod- 



uct of as described in Section 4.2, a change of basis can be accomplished 
using the appropriate unitary transformation: Hadamard transformation (de- 
noted H) changes between basis and basis, and Y = Rx{—'k/2) = 
exp(zc7^7r/4) transforms from basis to basis and from to . In 
matrix form, 



H 



1 



1 1 
1 -1 



Y 



1 

71 



1 i 
i 1 



Circuits of this form are the basis for the construction of the molecular unitary 



propagator as illustrated in Table Al where the circuit representations are 
given. 



5. The phase estimation algorithm 

In this section, we describe how to obtain molecular energies given the time 
evolution of the molecular Hamiltonian described above. The key idea is to 
Fourier transform the oscillating phase, (^(0)|V'(t)) = exp(— iSt), to obtain 
the electronic energy. The eigenenergy is converted into relative phases. The 
relative phase can then be measured using the quantum phase estimation 
algorithm (PEA). As the phase is measured the input state partially collapses 
to the set of states consistent with the measurements obtained up to that 
point. 

To determine an eigenvalue associated with an eigenstate, consider the phase 
of an eigenstate of the Hamiltonian, evolving dependent on a register qubit: 

|0)|V'n) +e-*^*|l)|V'„) = |0)|V^„) +e-*^"*|l)|V'„). (8) 

By letting En = 27r((/> — K)/t where < (/> < 1 and K is an unobservable 
integer, the unknown eigenvalue becomes encoded in the relative phase of the 
register qubit quantum state as |0) + g-^'^^^'^"-^) |1) [3 El [33]. Then, cf) is 



estimated using the phase estimation technique detailed in Subsection |5.1 

Once (j) has been estimated, it must be inverted to obtain the energy. Given 
bounds for the energy eigenvalue, [Emin, Emax), the time of evolution is se- 
lected as t = 2TT/{Emax — Emin) = 2tt / u and an energy shift of Eg is used to 
ensure K = {Eg — Emin)/^ is an integer. The energy shift, Eg, is effected by 
a gate on the register qubit which applies a phase to the qubit if it is in state 
|1) and does nothing otherwise. Using these parameters the measured value of 
(j) corresponds to the value of the energy, E^j, = ijj{4> — K) + Eg. 
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The PEA can be accompHshed in several different ways depending on the 
technology available and in the recent experimental implementations, the it- 
erative phase estimation algorithm [321 was used. We present it here 
for completeness. 



5.1. Iterative Phase Estimation 

The goal of the phase estimation algorithm is to obtain the best L bit estimate 
of a phase (p. The readout will be values of or 1 from the fractional binary 
expansion of (p: 



0.joiii2 • • • JL-I 



+ 



+ 



J2 



+ 



JL-I 

2^ 



(9) 



Each jfc is zero or one. The relative phase generated by the unitary prop- 
agator U{t) = exp{—iHt) applied as in ([s]) is defined as exp(27rii?!>). Per- 
forming U{2^t) results in a relative phase exp(27ri(2^0)). Observing that 
2^(j) = joji . . ■ jk-i-jkjk+i • • • we can finite binary expansion of length L de- 
terministically. 



First, 



IS 



implemented resulting in relative phase 



exp(27ri jo ■ ■ ■ Jl-i-Jl) = exp(27ri Jl/^). Since ji is or 1, the relative 
phase is either +1 or -1, respectively. The Hadamard transform described 
in Section 4.2.1 distinguishes |0) + |1) from |0) — |1) allowing a projective 
measurement to identify j^. 

To obtain Jl-i and more significant bits, the iterative method uses the gate 
Sk to deterministically obtaining jk - Since data is read from the least significant 
digit first, the counter-rotation Sk is a function of previously obtained bits. The 
form of the Sk gate is: 



1 




with = exp 



2Tri 



L-k+l . 

y-^ Jfc+z-i 

1=2 



2' 



This gate removes the L — k — 1 least significant digits so the state of quantum 
computer becomes (|0) + exp[— i7rjfc]|l))|^/;„) where jk is zero or one. Finally, 
effecting the Hadamard transformation leads to \jk)\ipn) and measurement of 
the register in the {|0), |1)} basis yields the value of jk- When the binary 
expansion of (j) is length L, the measurements are deterministic, otherwise the 
remainder causes some randomness but does not significantly affect the results. 
A complete error analysis can be found in Refs. |13| I35j. 

As noted in Ref. [37], in phase estimation algorithm the number of uses of 
U{to) for a fixed time Iq scales exponentially with the number of bits desired 
from (j). This is a consequence of the Fourier uncertainty principal; the more 
information required in the frequency domain the longer the propagation time. 
When the unitary is decomposed into gates, this means an exponential increase 
in the gates is required for an exponential increase in the precision of the 
measurement. 

To produce the state indicated in Q, the evolution of the eigenstate must be 
dependent on the register qubit requiring that the construction of the unitary 



evolution operator described in Section 4.2 be modified. The constructions 



listed in Table AT only need to be slightly modified; since the underlying 



constructions rely on gates, changing these rotations into controlled Rz 
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H 



H 



2k 



3k 



Figure 2. Iterative phase estimation circuit for obtaining the k^^ bit. The phase is represented using 
a binary decimal expansion (f) = 0.jojij2j3 ■ ■ ■ jn- To gain precision in the algorithm, one adjusts the 
parameterized gate Sk according to the values of all previous obtained bits. This (7 is a representation 
of the propagator of the system being simulated. Acting on the eigenstate with this operator advances 
the phase of the state allowing bit to be read out dcterministically. 



rotations (|1)(1| Rz + |0)(0| (8> 1) is sufficient to make the entire unitary 
dependent on the readout qubit. 

To obtain ground state energies, PEA reUes on the assumption that the 
input wave function has signification overlap with the ground state. Since each 
qubits represents the occupancy of molecular orbitals in the A^-electron wave 
function, the HF guess for the ground state requires no superposition 

of states and is thus straightforward to prepare. For \{'4'o^\'4'E^^)\ = 1 ~ 
where e is small, the phase estimation algorithm can be applied to retrieve an 
estimate of the ground state energy. Simultaneously, the state of the system 
will collapse to \'4^q^^) when measured in the H^^^ basis (via PEA) with high 
probability [32l[33]. If the Hartree-Fock guess is insufficient, more sophisticated 
state preparation procedures exist and these were reviewed in Ref. [12]. The 
adiabatic scheme for state preparation [7] is analyzed in the following section. 



6. Adiabatic state preparation 

When the initial state prepared as input has low overlap with the ground 
state, the initial state must be improved so that PEA will collapse the input 
to the correct state with greater efficiency. This section explains the method of 
preparing an input state to the simulation algorithm using adiabatic quantum 
computation [71 1381 139]. First, consider the Hartree-Fock wave function as an 
approximation to the FCI wave function. This is the output of an classical 
algorithm returning, in polynomial time, a computational basis state where 
the qubits which correspond to occupied molecular orbitals are in the state 
|1) with the remaining qubits in state |0). To increase overlap of the wave 
function, after the system is prepared in state \ipQ^), the Hamiltonian Hpci is 
slowly applied and the actual ground state is recovered by adiabatic evolution. 
Consider a smooth one-parameter family of adiabatic path Hamiltonians, 



H{s) = {1-s)Hhf + sHfci, 



(10) 



for monotonic s G [0, 1]. This was the adiabatic path originally proposed by 
us in Ref. [7|. Other paths may be used as in Ref. [ID] but in this study we 



restrict our attention to evolution of the form in Eq. (10). 



Let the instantaneous energies of H{s) be given by the sequence. 



Eois)<Ei{s) <■■■ <EN-iis), 



fill 



then the adiabatic state preparation procedure is efficient whenever the total 
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run time, T, satisfies the fohowing: 

r» min (El(s)-Eo(s))~^ (12) 

0<s<l 

according to known results relating the adiabatic theorem to complexity the- 
ory |38j . After adiabatic evolution, the state of the system is which 
is the ground state of the molecular Hamiltonian Hq'~'^ . Modified versions 
of this procedure exist using decoherence to achieve faster evolution and are 
discussed in Refs. [iHli2]. 

Assume that the adiabatic evolution induced transitions into higher energy 
states and so the un-normalized state of the system is + A|A;), where 

\k) ^Ti. = 1 — \4'o'~'^){'4'o^^\- While the error in the wave function is linear in 
A < 1, the overestimate of the energy in the expectation value {H^^^) is only 
quadratic. 

It is unclear how this method will scale. It is possible to prepare a desired 
ground state efficiently provided that the gap between the ground and excited 
states is sufficiently large [3S]. This depends on the adiabatic path taken. 
Finding the ground state energy of a random Hamiltonian, even for simple 
models, is known to be complete for the quantum analogue of the class NP |43[ 

E]. 

There are other ways to perform state preparation, for example, by going 
beyond Hartree-Fock theory |45l |l6] . A broader discussion of state preparation 
for quantum simulation can be found in Refs. [12l U?]. Recently, in Ref. [38] , 
the effects of initial states for the phase estimated quantum simulation CH2 
molecules were studied for a variety of geometries and eigenstates using initial 
guesses obtained via multi-configuration approaches. 



7. Simulating the hydrogen molecule 

To illustrate the algorithmic details of a scalable simulation of quantum sys- 
tems, the hydrogen molecule in a minimal basis is used as an instructive exam- 
ple. The minimal basis is the minimum number of spatial- functions needed to 
describe the system and in the case of H2, one spatial- function is needed per 
atom denoted ipHi and (pH2- In this simple case, the Hartree-Fock procedure 
is not necessary as the molecular spatial-orbitals are determined by symmetry 
and are given by ipg = ipui + ^H2 and (pu = ^hi — fH2- These two spatial 
functions correspond to four orbitals that will be identified as: 

\xi) = \^g)\a), \x2) = \^gm, \x3) = \^u)\a), \x4) = \y^um. (13) 

The form of the spatial function is determined by the basis set used. The 
ST0-3G basis is a common minimal basis that approximates a single electron 
spatial Slater type orbitals (STO), with a contraction of three real Gaussian 
functions [1]. Using this orbital basis, the spatial integrals of the Hamiltonian 
were evaluated in Table [2] for H2 at bond distance 1.401000 atomic units 
(7.414-10-11 m). 

Considering H from Eq. ([T]) as H^^^ + H^"^^ we have, 

7/(1) = hiia[ai + /i22a2"2 + ^330303 + h^^ala^ (14) 
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Table 2. The one-electron and two-electron integrals defined in 
Eqs. |3j and ijl} arc evaluated using the molecular spatial orbitals 
obtained from a restricted Hartrce-Fock calculation at an internu- 
clear distance of 1.401000 atomic units (7.414-10"^^ m) [iol . 



Integrals 


Value (a.u.) 


/ill = /l22 
/133 = /144 


-1.252477 
-0.475934 


/11221 = '12112 


0.674493 


'13443 = '14334 


0.697397 


'11331 = '1144I = '12332 = '14224 
= '13I13 = '13223 = '14114 = '14224 

'11313 = '12424 = '13241 = '11423 = 'll243 


0.663472 
0.181287 



The following circuit applies the single-electron propagator for a time t: 



X2 
X3 
X4 







< 




^T(hiit)^ 








1 1 1 






1 1 K'^33l) 1 





T{h33t) 



The gate T is defined as: 



j{e) 



(15) 



The two electron Hamiltonian can also be expanded. As electrons are indis- 
tinguishable, 



hpqrs — 



the two electron Hamiltonian can 



Xp(x2)Xg(xi)Xr(xi)Xs(x2) 



and apa^qa^as = aJapOsaj. 



ri2 

be simplified as: 

H'^'^^ = /li22iaia2'^2«l + ^3443«3a4a4a3 + /ll44ia| 040401 -I- /i233202a30302 

+ (/11331 - /ii3i3) o| 0^^0301 + (/i2442 - /12424) a\a\aia2 

-f-3f?(/il423)(«l 040203 + 0^^0^0401) -I- K(/li243) (01020403 -I- 0^^0^0201) 
-|-9(/ll423)(o| 0^0203 - 0^0^0401) -1- 9(/ll243) (010^0403 - 0^50^020^6) 

The first six terms, 



''qpsr 1 



/ll22lo| 030201 + /13443O3O4O4O3 -I- /li44io| 040401 



+ /l2332a2030302 + (/ll331 " ^^ISls) o{ 0^0301 + (/l2442 " /12424) 0^0^0402 



can be simulated using the system Hamiltonian that employs only commut- 
ing two-local terms described in Section |4j Notice after the Jordan- Wigner 
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transform of the relevant operator we have: 



13 



p<q 



V p<q 



pqqp "-pqpq 



hpnpqSa„(7„) 1 



■ ^ E I y^X^pqqp hpgpgda^aj j CT^ 

1 \pf^q J 



p<q 



{hpqqp hpqpqSa^a,) zz 



aXq. (17) 



The factor of 1 /2 is accounted for because the indistinguishably of the electrons 
reduces the summatfon in Eq. ([T]). 

Following Eq. ([l7]), let 6 = (1/4) X]p<q(^pggp " hpqpq6a^aj and 9p = 
Yliq-pj!zqi^pqqp ~ hpqpq^(Tj,aJ- Then following circuit illustrates the one and two- 
local interactions required to implement Eq. (17): 



- jT(et)| 1 i- 



PEA 

XI 1 -Rz(flit) [ 

X2 

X3 

X4 



fl,(fl4t) 



Defining ijpg 
picted as: 



^i^T-pqqp ~ hpqpg6a-paj, the three local interactions can be de- 



^ "z('?23*T] 'l) q) I Hz (US') [ 4- 



-$- | Hz(l34'0] (P CP ^^"^(124') | q) CP I flz('714') \ A 



Each term commutes thus can be realized in any order. 
The remaining terms are strictly real leaving: 



hi423{o\a\a2a3 + 03030401) + (/ii243)(aia2"4a3 + O3O4O2O1). 

Since the orbitals are real, the integrals /ii243 = /11423 are equivalent. Therefore, 
we are left with the task of simulating 2/11423(0! O4O2O3 + 03030401). 

Consider the general term OpOgOfOs+oJoJogOp. Due to the anti-commutation 
rules, all sets of operators corresponding to a set of four distinct spin-orbitals, 
{p, q, r, s), are simulated using the same circuit. This is due to the fact that 
the Jordan- Wigner of each operator generates same set of operators (namely, 
the eight combinations involving an even number of and operators). 
However, depending on if or a~ is used each term of spin operators will 
have a different sign. If we define: 



^ — {hpgrs^apCTs^cTqa,, ^qprs^crpa^^cTqas) 
^ — (hpsqr^apCTr^(Jq(Ts ^spqr^cTpUq^arCTs^ 
— {hprsq^aj,CFq^(Tr(Js f^prqs^(Tpas^(Tqar) 



(18) 
(19) 
(20) 
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(a) Plot of gates to simulate the H2 (b) Plot of relative error of approxima- 
versus time step used in the first order tion as a function of gates used. 
Trotter approximation. 

Figure 3. The unitary propagator corresponding to this Hamiltonian is approxi- 
mated using a first order Trotter decomposition and these graphs provide analysis 
of the Trotter error and the number of gates used at each Trotter number, T„. The 
unitary propagator is simulated by applying each small term of the second quantized 
Hamiltonian for small time steps dt and repeating the process Tn = t/dt times. As 
dt decreases the error in the approximation decreases at the expense of more gates. 
Zero error represents the eigenvalue of the Hamiltonian of H2 in the minimal basis at 
a separation of 1.4 atomic units. The horizontal line of (b) represents the threshold 
for energy error of 10"'* atomic units. 



then 



s-1 



(21) 



,k=p+l k=r+l 



Applying this to the hydrogen molecule, observe that h^^^ = —h^"^^ and h^^^ = 
indicating that only the terms {(j^cr^cr^cr^, cr^cr^fj^cr^, (T^cr^a^'cr^, cj'^cr^cj^fT^} 
must be considered. The resulting quantum circuit is illustrated in Table |A3[ 

To assess the Trotter error, we simulated this circuit using a classical com- 
puter using the first-order Trotter decomposition. The pseudo-code for the H2 
simulation is given in Appendix |A] and the results are summarized in Fig. 3. 
Although the gates increase with the Trotter number, reducing the Trotter 
error of the dynamics is only practical if the measurement is precise enough 
to detect such errors. Thus, in practice, there is a balance between the Trotter 
number selected and the number of bits to be obtained by the measurement 
procedure. Finally, as noted in section 4.1 higher-order Trotter decomposi- 



tions also allow for more precise simulation with larger time steps. In fact, 
given the largest eigenvalue of the simulated system as A the number of gates 
for simulation time t is 0{tX) using the appropriate order of Trotter decom- 
position [28]. 



8. Conclusions 



In this paper, we mapped the full configuration interaction (FCI) method 
from quantum chemistry into a quantum algorithm. We reviewed the elec- 
tronic structure problem, techniques of creating the simulated propagator. 
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and explicitly illustrated this construction for H2 for a single time step of a 
first-order Trotter expansion. 

Applicability of quantum simulation comes down to the ability to propagate 
the simulated system with a specified error tolerance. Since phase estimation 
is essentially a Fourier transform of the frequency of phase oscillations (which 
are proportional to the eigenenergy) to obtain more precise determination of 
the frequency, a longer propagation time is necessary. The longer simulation 
time requires more manipulations of the computational system. 

In the coming future, small scale experiments such as the simulation of 
the circuits we have presented for H2 on a quantum computer will likely be 
possible. Experimental realizations of quantum chemistry on quantum devices 
have only recently been achieved [T71 [18] . We hope that the present paper will 
continue the interest by giving an example of a scalable construction of the 
unitary propagator for the H2 molecule in an explicit form, which poses the 
next logical challenge for experimental realization. 
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Appendix A. Example program for a quantum computer to calculate 
electronic energies for the hydrogen molecule 

The time evolution of input states is measured to extract eigenvalues of the 
system. A state preparation procedure must be used if a particular state is 
desired (e.g. the ground state). Time evolution is implemented in such a way 
that a measurable relative phase is generated on a register qubit. This register 
qubit is then measured to obtain a bit of information about the phase. The time 
evolution is then effected for half as long and again the register is measured 
to obtain the next digit. This process is repeated until all digits are read out 
and then the phase can be converted to an energy eigenvalue following the 
prescription of Section [5] 

In this appendix, we explicitly spell out the gates that would be necessary 
for effecting a single time step of the evolution generated by the Hydrogen 
molecule's Hamiltonian. The time step idea is to alternatively apply each of 
the non-commuting unitary operators for a fraction of the total time and 
repeat the time steps until the total time has been reached. As the length 
of the time steps goes to zero the error of the approximation vanishes. For 
constructing long time evolution the time steps should be composed following 



the Trotter prescription 4.1 



The decomposition of the propagator of H2 Hamiltonian is given using 
standard quantum gates. The entire unitary operator is applied if the reg- 
ister qubit is in state Qubits are named: Register, Ql, Q2, Q3, Q4. 
The integrals hij and hijki are parametrized by the nuclear configuration 
and are given in Table [2] for the for the equilibrium geometry. Addition- 
ally, we define riij = {hijji — hjiji)/A. Note only /11313 and /i2424 are the 
only non-zero terms due to spin orthogonality. Following the main text, let 

Q = (1/4) '^p^q{hpqqp — hpqpq5a^cr^) and Op = q-p-Lq{hpqqp — hpqpqdcr^ag) ■ 
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Gate 


Target qubit 


Control qubit 


Parameter 


Hadamard 


Register 






Single Electron Operators 






cPhase 


Ql 


Register 


hut 


cPhase 


Q2 


Register 


h22t 


cPhase 


Q3 


Register 




cPhase 


Q4 


Register 


/l.44t 


Two Electron Operators: number- 


■NUMBER operator 




Phase 


Register 




Qt 


cR^ 


Ql 


Register 


-Oit 


cR^ 


Q2 


Register 


-e2t 


cRz 


Q3 


Register 


-03t 


cRz 


Q4 


Register 




cNot 


Q4 


Q3 




cR^ 


Q4 


Register 




CiNot 


Q4 


Q3 




cNot 


Q4 


Q2 




cRz 


Q4 


Register 


2n2At 


cNot 


Q4 


Q2 




cNot 


Q4 


Ql 




cRz 


Q4 


Register 


2niit 


cNot 


Q4 


Ql 




cNot 


Q3 


Q2 




cRz 


Q3 


Register 


2n23t 


cNot 


Q3 


Q2 




cNot 


Q3 


Ql 




cRz 


Q3 


Register 


2ni3i 


cNot 


Q3 


Ql 




cNot 


Q2 


Ql 




cR^ 


Q2 


Register 


2ni2t 


cNot 


Q2 


Ql 




Two Electron Operators: excitation- excitation operator 


XXYY 








Hadamard 


Ql 






Hadamard 


Q2 






Rx 


Q3 




-tt/2 


Rx 


Q4 




-7r/2 


cNot 


Q2 


Ql 




cNot 


Q3 


Q2 




cNot 


Q4 


Q3 




cR^ 


Q4 


Register 


-t{hi423 + hi243)/4: 


cNot 


Q4 


Q3 




cNot 


Q3 


Q2 




cNot 


Q2 


Ql 




Rx 


Q4 




7r/2 


Rx 


Q3 




7r/2 


Hadamard 


Q2 







continued on next page 
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Gate 


TarP'P'l" miniti 

-L OjJ- (j VJ LI Ul \j 


(innf.rol nnhif 


Pa r Q rn pfPT' 


Hadamard 


Ql 






YYXX 










Ql 




-tt/2 


R-x 


Q2 




-tt/2 


Hadamard 


Q3 






Hadamard 


Q4 






cNot 


Q2 


Ql 




cNot 


Q3 


Q2 




cNot 


Q4 


Q3 




cR^ 


Q4 


Register 


-t{hu23 + /il243)/4 


cNot 


Q4 


Q3 




cNot 


Q3 


Q2 




cNot 


Q2 


Ql 




Hadamard 


Q4 






Hadamard 


Q3 








Q2 




7r/2 


Rx- 


Ql 




7r/2 


XYYX 








Hadamard 


Ql 






Rx 


Q2 




-tt/2 


Rx 


Q3 




-7r/2 


Hadamard 


Q4 






cNot 


Q2 


Ql 




cNot 


Q3 


Q2 




cNot 


Q4 


Q3 




cR^ 


Q4 


Register 


i(/i-1423 + /ll243)/4 


cNot 


Q4 


Q3 




cNot 


Q3 


Q2 




cNot 


Q2 


Ql 




Hadamard 


Q4 






Hadamard 


Q3 






Rx 


Q2 




7r/2 


R-x 


Ql 




tt/2 


YXXY 








Rx 


Ql 




tt/2 


Hadamard 


Q2 






Hadamard 


Q3 






Rx 


Q4 




tt/2 


cNot 


Q2 


Ql 




cNot 


Q3 


Q2 




cNot 


Q4 


Q3 




cR^ 


Q4 


Register 


t{hi423 + /il243)/4 


cNot 


Q4 


Q3 




cNot 


Q3 


Q2 




cNot 


Q2 


Ql 




Rx 


Ql 




-tt/2 


Hadamard 


Q2 






Hadamard 


Q3 








Ql 




-n/2 
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Table Al. The quantum circuits corresponding to evolution of the listed Hermitian second-quantized operators. Here, p, q, 
r, and s arc orbital indices corresponding to qubits such that qubit state |1) indicates an occupied orbital and |0) indicates 
unoccupied. It is assumed that the orbital indices satisfy p > g > r > s. These circuits were found by performing the Jordan- 
Wigner transformation given in Eqs. jsbj and (jSaJ and then propagating the obtained Pauli spin variables |30| . In each 
circuit, 6 — 0{h) where h is the integral preceding the operator. Gate T{0) is defined by T|0) — |0) and T|l) — exp{ — 20)|1), 
G is the global phase gate given by exp( — and the change-of-basis gate Y is defined as Rx{—7r/2). Gate H refers to the 
Hadamard gate. For the number-excitation operator, both M — Y and M — H must be implemented in succession. Similarly, 
for the double excitation operator each of the S quadruplets must be implemented in succession. The global phase gate must 
be included due to the phase-estimation procedure. Phase estimation requires controlled versions of these operators which 
can be accomplished by changing all gates with ^-dependence into controlled gates dependent on register qubits. 



Second quantized operators 



Circuit 



Number 
operator 



hpp dp dp 



T{0) 



Excitation 
operator 



hp,(a+a<j+a+ap) 



-i — L 



Coulomb and 
exchange operators 



Number-excitation" 
operator 



i- 



G(e) - Rz{e) 



Rz(0) 



Rz(0) -e- 



p -Tm}-^ 



q + i 
q 

9-1 



I where M={H,Y} 















( 


^ — 





) 



] — ^ 




Rz{d) 




^ [ 











Double excitation 
operator 



(at a^ds 



Ml 



M2 



Ma 



i- 


1 














— ( 


>^ 


) — 




1 1 





m{ 



M^ 



Ml 



where (Mi,M2,M3,M4)= 
{{H,H,H,H), (Y,Y,Y,Y), 
(H,Y,H,Y),(Y,H,Y,H), 
(Y,Y,H,H),(H,H,Y,Y), 
(Y,H,H,Y), (H,Y,Y,H)} 
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Notation: 



^The spin variable representation of this operator depends on whether q lies in the range p to r or outside of it. 



December 21, 2010 



1:24 



arXiv: 1001.3855 



h2.arxiv2 



arXiv:1001.3855 



21 



C5 _>> 

^ fvJ CD . 



s .s 



r-( O tn ^ O 

Cl, g u cfl o 

j3 bO c3 «y 

CO +J [Z] S-I ^ 



C ^ --t^ 



0) o 
(ft pi 



■2 ^ I g 



4j T3 



3 S' 



S oi a 

+j 

Ml S " ? 

H 52 c „ 

■g ^^ QJ 

F o c ^ 



9- e 



0) 03 



*2 M M 



• P c B H 



tt -a 



:li 

B . 
.3 



OS :j *3 B 

H ? » I 

O a a, S 3 

b ^ QJ O 

tS B ° " O 



a. 2 



'is? 



CD +2 

£ « 5 P 

< .S a - 

0, -Q (ft 2 

H g 8 .2 



J3 



O 
u 

CD 



CD 
P 




e 

H— 0< 

e 

-I— s. 

e 
a. 



CS 



a. 



C3 
a. 



O 

CD 

3 S 
^ a 
^ O 



.9 S 
H O 



CD 

a 
O 

a 
o 



CD 

a 



O O 



CD 
I 



CD 



.2 S 

c^ cS 



CD 

X 

H 



o 

P H O 



3 oj 



cs 



B 



o 



.2 O 



A 



December 21, 2010 1:24 arXiv: 1001.3855 h2.arxiv2 



22 Whitfield et. al 



a III ° 



a "u 
■5 c 



g 5 % 



O tD 



k| £ 

KIN C3 C! 

O O 

CD W _ 

I- oj d 

Si „ 

^ O u 

O O 0) 

?3 OJ 

S £ s 
s ^ 

c S 1 



o - " 

« a 

" CD 

^ _aj 



* I 

0) 0) - 

1g £ ^ 

.2 *^ j:j 

u t+H tn 

0-- 

m K 

m jc CIS 

* "S > 



CO 



4j CD o 

S eg 



bo 



H B £ 

J5 c o 



S H 

a c 

s s " 

J3 1^ 

^ -0 

^ ^ N 

■5 > 

u t- . 

O T) II 



gild 
<■ .ti + ^ 

! B 
H H M 



m @ @ m 




-9 

a; 




a; 



i @ 5 5 



a, 



